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We present the computation of the full system of Boltzmann equations at second-order describing 
>P ■ the evolution of the photon, baryon and cold dark matter fluids. These equations allow to follow 

the time evolution of the Cosmic Microwave Background (CMB) anisotropies at second-order at all 
angular scales from the early epoch, when the cosmological perturbations were generated, to the 
present through the recombination era. This paper sets the stage for the computation of the full 
^~^' second-order radiation transfer function at all scales and for a a generic set of initial conditions 

C^ \ specifying the level of primordial non-Gaussianity. In a companion paper, we will present the 

computation of the three-point correlation function at recombination which is so relevant for the 
issue of non-Gaussianity in the CMB anisotropies. 



o 
o 

(N 



in 
> 



6 

> 

X 



PACS numbers: 98.80.Cq DFPD 06/A/05 



■^ ■ I. INTRODUCTION 

\^ ' Cosmological inflation [1] has become the dominant paradigm to understand the initial conditions for the CMB 
f^ , anisotropies and structure formation. In the inflationary picture, the primordial cosmological perturbations are 
created from quantum fluctuations "redshifted" out of the horizon during an early period of accelerated expansion of 
the universe, where they remain "frozen". They are observable as temperature anisotropies in the CMB. This picture 
has recently received further spectacular conflrmations by the Wilkinson Microwave Anisotropy Probe (WMAP) three 
O ' year set of data [2]. Since the observed cosmological perturbations are of the order of 10~^, one might think that 
flrst-order perturbation theory will be adequate for all comparisons with observations. This might not be the case, 
^ ' though. Present [2] and future experiments [3] may be sensitive to the non-linearities of the cosmological perturbations 
1^ • , at the level of second- or higher-order perturbation theory. The detection of these non-linearities through the non- 
Gaussianity (NG) in the CMB [4] has become one of the primary experimental targets. 

There is one fundamental reason why a positive detection of NG is so relevant: it might help in discriminating among 
""jh ' the various mechanisms for the generation of the cosmological perturbations. Indeed, various models of inflation, firmly 
Cu ^ rooted in modern particle physics theory, predict a significant amount of primordial NG generated either during or 
immediately after inflation when the comoving curvature perturbation becomes constant on super- horizon scales [4]. 
While single-fleld [5] and two(multi)-fleld [20] models of inflation predict a tiny level of NG, 'curvaton-type models', 
in which a signiflcant contribution to the curvature perturbation is generated after the end of slow-roll inflation by the 
perturbation in a field which has a negligible effect on infiation, may predict a high level of NG [7] . Alternatives to the 
curvaton model are those models characterized by the curvature perturbation being generated by an inhomogeneity 
in the decay rate [8, 9], the mass [10] or the interaction rate [11] of the particles responsible for the reheating after 
infiation. In that case the reheating can be the first one (caused by the scalar field(s) responsible for the energy density 
during infiation) or alternatively the particle species causing the reheating can be a fermion [12]. Other opportunities 
for generating the curvature perturbation occur at the end of inflation [13], during preheating [14], and at a phase 
transition producing cosmic strings [15]. 

Statistics like the bispectrum and the trispectrum of the CMB can then be used to assess the level of primordial NG 
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on various cosmological scales and to discriminate it from the one induced by secondary anisotropics and systematic 
effects [4, 16, 17, 18]. A positive detection of a primordial NG in the CMB at some level might therefore confirm 
and/or rule out a whole class of mechanisms by which the cosmological perturbations have been generated. 

Despite the importance of NG in the CMB anisotropics, little effort has ben made so far to provide accurate 
theoretical predictions of it. On the contrary, the vast majority of the literature has been devoted to the computation 
of the bispectrum of cither the comovig curvature perturbation or the gravitational potential on large scales within 
given inflationary models. These, however, are not the physical quantities which are observed. One should instead 
provide a full prediction for the second-order radiation transfer function. A firt step towards this goal has been 
taken in Ref. [19] (see also [20]) where the full second-order radiation transfer function for the CMB anisotropics 
on large angular scales in a flat universe filled with matter and cosmological constant was computed, including the 
second-order generalization of the Sachs- Wolfe effect, both the early and late Integrated Sachs- Wolfe (ISW) effects 
and the contribution of the second-order tensor modes. 

There are many sources of NG in the CMB anisotropics, beyond the primordial one. The most relevant sources are 
the so-called secondary anisotropies, which arise after the last scattering epoch. These anisotropics can be divided 
into two categories: scattering secondaries, when the CMB photons scatter with electrons along the line of sight, 
and gravitational secondaries when effects are mediated by gravity [21]. Among the scattering secondaries we may 
list the thermal Sunyaev-Zeldovich effect, where hot electrons in clusters transfer energy to the CMB photons, the 
kinetic Sunyaev-Zeldovich effect produced by the bulk motion of the electrons in clusters, the Ostrikcr-Vishniac effect, 
produced by bulk motions modulated by linear density perturbations, and effects due to reionization processes. The 
scattering secondaries are most significant on small angular scales as density inhomogeneitics, bulk and thermal 
motions grow and become sizeable on small length scales when structure formation proceeds. 

Gravitational secondaries arise from the change in energy of photons when the gravitational potential is time- 
dependent, the ISW effect, and gravitational lensing. At late times, when the Universe becomes dominated by the 
dark energy, the gravitational potential on linear scales starts to decay, causing the ISW effect mainly on large angular 
scales. Other secondaries that result from a time dependent potential are the Rees-Sciama effect, produced during 
the matter-dominated epoch at second-order and by the time evolution of the potential on non-linear scales. 

The fact that the potential never grows appreciably means that most second order effects created by gravitational 
secondaries are generically small compared to those created by scattering ones. However, when a photon propagates 
from the last scattering to us, their path may be deflected because of the gravitational lensing. This effect does not 
create anisotropies, but it only modifies existing ones. Since photons with large wavenumbers k are lensed over many 
regions {^ k/H, where H is the Hubble rate) along the line of sight, the corresponding second-order effect may be 
sizeable. The three-point function arising from the correlation of the gravitational lensing effect and the ISW effect 
generated by the matter distribution along the line of sight [22, 23] and the Sunyaev-Zeldovich effect [24] are large 
and detectable by Planck [25]. 

Another relevant source of NG comes from the physics operating at the recombination. A naive estimate leads to 
think that these non-linearities are tiny being suppressed by an extra power of the gravitational potential. However, 
the dynamics at recombination is quite involved because all the non-linearities in the evolution of the baryon-photon 
fluid at recombination and the ones coming from general relativity should be accounted for. This complicated dynamics 
might lead to unexpected suppressions or enhancements of the NG at recombination. A step towards the evaluation 
of the three-point correlation function has been taken on Ref. [26] where it was computed in the in so-called squeezed 
triangle limit, when one mode has a wavelength much larger than the other two and is outside the horizon. 

In this paper wc present the computation of the full system of Boltzmann equations at second-order describing the 
evolution of the photon, baryon and CDM fluids, neglecting polarization. These equations allow to follow the time 
evolution of the CMB anisotropies at second-order at all angular scales from the early epoch, when the cosmological 
perturbations were generated, to the present through the recombination era. This paper sets therefore the stage for the 
computation of the full second-order radiation transfer function at all scales and for a a generic set of initial conditions 
specifying the level of primordial non-Gaussianity. In a companion paper [27], we will present the computation of the 
three-point correlation function at recombination which is so relevant for the issue of NG in the CMB anisotropies. 
Of course on small angular scales, fully non-linear calculations of specific effects like Sunyaev-Zel'dovich, gravitational 
lensing, etc. would provide a more accurate estimate of the resulting CMB anisotropy, however, as long as the leading 
contribution to second-order statistics like the bispectrum is concerned, second-order perturbation theory suffices. 

The paper is organized as follows. In Section II we provide the second-order metric and corresponding Einstein 
equations. In Section HI the left-hand-side of the Boltzmann equation for the photon distribution function is derived 
at second-order. The collision term is computed in Section IV. In Section V, we present the second-order Boltzmann 
equation for the photon brightness function, its formal solution with the method of the integration along the line of 
sight and the corresponding hierarchy equations for the multipole moments. Sections VI contains the derivation of the 
Boltzmann equations at second-order for baryons and Cold Dark Matter (CDM). Section VII contains the expressions 
for the energy momentum tensors and, finally. Section VIII contains our summary. 



In performing the computation presented in this paper, we have mainly followed Ref. [28] (in particular chapter 
4) where an excellent derivation of the Boltzmann equations for the baryon-photon fluid at first-order is given. Since 
the derivation at second-order is straightforward, but lenghty, the reader might benefit from reading the appropriate 
sections of [28]. At the end of the paper we have provided a Table which summarizes the many symbols adopted 
appearing throughout the paper. 

II. PERTURBING GRAVITY 

Before tackling the problem of interest - the computation of the Boltzmann equations for the baryon-photon and 
CDM fluids - we first provide the necessary tools to deal with perturbed gravity, giving the expressions for the Einstein 
tensor perturbed up to second-order around a flat Friedmann- Robertson- Walker background, and the relevant Einstein 
equations. In the following we will adopt the Poisson gauge which eliminates one scalar degree of freedom from the 
goi component of the metric and one scalar and two vector degrees of freedom from g^ . We will use a metric of the 
form 

ds"^ = a^ir]) [-e'^'^dr]'^ + 2uj^dx'dr] + (e^^^^ + Xr])dx'dx^] , (2.1) 

where 0(77) is the scale factor as a function of the conformal time 77, and uji and Xij ^re the vector and tensor 
peturbation modes respectively. Each metric perturbation can be expanded into a linear (first-order) and a second- 
order part, as for example, the gravitational potential $ = <I>'^^ + $'^^^2. However in the metric (2.1) the choice of 
the exponentials greatly helps in computing the relevant expressions, and thus we will always keep them where it 
is convenient. From Eq. (2.1) one recovers at linear order the well-known longitudinal gauge while at second-order, 
one finds $(2) ^ ^(2) _ 2(0(i))2 and *(2) ^ ^(2) _^ 2(V'(i))2 where 0(1), V^*^^ and 0(2), t^i^) (with 0(i) = $(1) and 
■0(1) = vl>(i)) are the first and second-order gravitational potentials in the longitudinal (Poisson) gauge adopted in 
Refs. [4, 29] as far as scalar perturbations are concerned. For the vector and tensor perturbations, we will neglect linear 
vector modes since they are not produced in standard mechanisms for the generation of cosmological perturbations (as 
inflation), and we also neglect tensor modes at linear order, since they give a negligible contribution to second-order 
perturbations. Therefore we take LVi and Xij to be second-order vector and tensor perturbations of the metric. 

Let us now give our definitions for the connection coefficients and their expressions for the metric (2.1). The 
number of spatial dimensions is n = 3. Greek indices (a,/3, ..., /i, j^, ....) run from to 3, while latin indices 
(a, 5, ..., i,j, fc, ....TO, n...) run from 1 to 3. The total spacetime metric g^i, has signature (— , -I-, +, -I-). The connection 
coefficients are defined as 



The Riemann tensor is defined as 
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The Ricci tensor is a contraction of the Riemann tensor 

J^pv = n ^^^, , (2.4) 

and in terms of the connection coefficient it is given by 

n> _ a pa _ a pa I pa p<T _ pa p<T /r, r\ 
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The Ricci scalar is given by contracting the Ricci tensor 

R = R^^. (2.6) 

The Einstein tensor is defined as 

Gp^ ~ Rf^i, — -g^j^ijR. (2.7) 



The connection coefficients 



For the connection coefRcients we find 
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B. Einstein equations 

The Einstein equations are written as G^i, — k^T^^, so that k^ — SirG-^, where Gn is the usual Newtonian gravitational 
constant. They read 

-2$ 

G^o = -^ [3H^ - 6H^' + 3(*')' - e'*+'* (9^*9^* - 2V2*)] = K^yO^ ^ (2.9) 

G\ = 2^ [a**' + (W - *') a*$] - -^V^Lu' + 4^2 _ 2_ _ = k2j,^^ ^ (2.10) 
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- 5 (9H + 5,-0 - ^ {d^ + d^u.') + ^ (Wx}' + ^Xf - iv^x}) = -'T) . (2.11) 

Here T^^ is the energy momentum tensor accounting for different components, photons, baryons, dark matter. We 
will give the expressions later for each component in terms of the distribution functions. 

III. THE COLLISIONLESS BOLTZMANN EQUATION FOR PHOTONS 

We are now interested in the anisotropics in the cosmic distribution of photons and inhomogeneities in the matter. 
Photons are affected by gravity and by Compton scattering with free electrons. The latter are tightly coupled to 
protons. Both are, of course, affected by gravity. The metric which determines the gravitational forces is influenced 
by all these components plus CDM (and neutrinos). Our plan is to write down Boltzmann equations for the phase- 
space distributions of each species in the Universe. 

The phase-space distribution of particles g{x^,P^,'q) is a function of spatial coordinates x', conformal time 77, 
and momentum of the particle P'' = dx^/dX where A parametrizes the particle path. Through the constraint 
p2 = gfii^P^P" = —m?, where m is the mass of the particle one can eliminate P" and g(a;% P^, 77) gives the number 
of particles in a differential volume dx^dx^dx^dP^dP^dP^ in phase space. In the following we will indicate the 
distribution function for photons with /. 

The photons' distribution evolves according to the Boltzmann equation 

f=aG[/], (3.1) 

drj 

where the collision term is due to scattering of photons off free electrons. In the following we will derive the left-hand 
side of Eq. (3.1) while in the next section we will compute the collision term. 



For photons wc can impose P^ = g^iyP'^P'^ = and using the metric (2.1) in the conformal time 77 we find 
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where we define 



From the constraint (3.2) 



p2 ^ g^^P^pl . (3.3) 
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P° = e-*f^ + 2c^,P"PM . (3.4) 

Notice that we immediately recover the usual zero and first-order relations P° = p/a and P" = p{l — $'^^))/a. 

The components P* are proportional to pn^, where n* is a unit vector with n^rii = SijU^n^ — 1. We can write 
P' = Cn^ , where C is determined by 

so that 

pi = ^„^ (e-2* + xkmn^ff^y"^ = ^n^e* f 1 - KkmU^'nA , (3.6) 

a a \ 2 J 

where the last equality holds up to second order in the perturbations. Again we recover the zero and first-order 
relations P* = pn^ /a and P* = pn'^{l + 'i>^^')/a respectively. Thus up to second order we can write 

P'> = e-'^^{1 + Lu,n') . (3.7) 
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Eq. (3.6) and (3.7) allow us to replace P° and P* in terms of the variables p and n*. Therefore, as it is standard in 
the literature, from now on we will consider the phase-space distribution / as a function of the momentum p = pn^ 
with magnitude p and angular direction n*, f = f{x^,p,n^,r]). 

Thus in terms of these variables the total time derivative of the distribution function reads 



df _ df df dx' df dp df dn' 
drj drj dx"^ drj dp drj dn^ drj 



(3.8) 



In the following wc will compute dx^ /drj, dp /drj and dn^ /drj: 
a) dx^/drj: 



From 



and from Eq. (3.6) and (3.7) 
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^ = n^e*+* (1 - oj, n-'" - ^Xk^n^nA . (3.10) 

b) dp/drj: 

For dp/drj we make use of the time component of the geodesic equation dP^ /dX = — F" oP^P'' where we can 
replace the derivative d/dX with (drj/dX) d/dt] = P" d/dr] so that 
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= -C^5?^> (3-11) 



Using the metric (2.1) we find 
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On the other hand the expression (3.7) of P*^ in terms of p and n* gives 
Thus Eq. (3.11) aUows us express dp/drj as 
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where in Eq. (3.12) we have replaced P° and P* by Eqs. (3.7) and (3.6). Notice that in order to obtain Eq. (3.14) we 
have used the following expressions for the total time derivatives of the metric perturbations 
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where we have taken into account that uji is already a second-order perturbation so that wc can neglect dn^ /drj which 
is at least a first order quantity, and we can take the zero-order expression in Eq. (3.10), dx"^ /drj = n^. In fact there 
is also an alternative expression for dp/drj which turns out to be useful later and which can be obtained by applying 
once more Eq. (3.15) 



1 dp ^ , d^ ^ , ^ . , , 1 / , 1 

p drj drj 2 ■' 



(3.17) 



c) dn^ /drj: 



We can proceed in a similar way to compute dn^/drj. Notice that since in Eq. (3.8) it multiplies df/dri^ 
which is first order, we need only the first order perturbation of dn^/drj. We use the spatial components of the 
geodesic equations dP^/dX = —T^^i^P'^Pl^ written as 



dP' 
drj 

For the right-hand side we find up to second-order 
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while the expression (3.6) of P* in terms of our variables p and n* in the left-hand side of Eq. (3.18) brings 
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Thus, using the expression (3.6) for P' and (3.4) for P'^ in Eq. (3.19), together with the previous result (3.14), the 
geodesic equation (3.18) gives the following expression dn^/drj (valid up to first order) 

^ = ($.fe + *,fc)n'=n^-$'*-*'\ (3.21) 

dry 

To proceed further we now expand the distribution function for photons around the zero-order value f^^' which is 
that of a Bose-Einstein distribution 

f^"Hp,v)^2 r^^ , (3.22) 

where T(ri) is the average (zero-order) temperature and the factor 2 comes from the spin degrees of photons. The 
perturbed distribution of photons will depend also on x* and on the propagation direction n* so as to account for 
inhomogeneities and anisotropics 

/(a;^p,n^r7) = /(")(p,,?) + /(l)(a;^p,n^ry) + i/(2)(x^p,n^,7), (3.23) 

where we split the perturbation of the distribution function into a first and a second-order part. The Boltzniann 
equation up to second-order can be written in a straightforward way by recalling that the total time derivative of a 
given i-th perturbation, as e.g. df^^> /drj is at least a quantity of the i-th order. Thus it is easy to realize, looking at 
Eq. (3.8), that the left-hand side of Boltzniann equation can be written up to second-order as 

df d/(^) , 1 d/(^) ^a/(o) d f^,,, ^ 1 ^(,A ^ dj^d_ / (,) ^ ^(,) ^ 1^(,) ^ 1^(,) 
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where for simplicity in Eq. (3.24) we have already used the background Boltzmann equation {df /drj)Y^> — 0. In 
Eq. (3.24) there are all the terms which will give rise to the integrated Sachs- Wolfe effects (corresponding to the terms 
which explicitly depend on the gravitational perturbations), while other effects, such as the gravitational lensing, are 
still hidden in the (second-order part) of the first term. In fact in order to obtain Eq. (3.24) we just need for the time 
being to know the expression for dp /drj, Eq. (3.17). 

IV. COLLISION TERM 



A. The Collision Integral 

In this section we focus on the collision term due to Compton scattering 

e(q)7(p) — e(q')7(p') , (4-1) 

where we have indicated the momentum of the photons and electrons involved in the collisions. The collision term 
will be important for small scale anisotropics and spectral distortions. The important point to compute the collision 
term is that for the epoch of interest very little energy is transferred. Therefore one can proceed by expanding the 
right hand side of Eq. (3.1) both in the small perturbation, Eq. (3.23), and in the small energy transfer. Part of the 
computation up to second-order has been done in Refs. [30, 31, 32] (see also [33]). In particular Refs. [30, 31] are 
focused on the effects of reionization on the CMB anisotropics thus keeping in the collision term those contributions 
relevant for the small scale effects due to reionization and neglecting the effects of the metric perturbations on the 
left-hand side of Eq. (3.1). We will mainly follow the formalism of Ref. [31] and we will keep all the terms arising 
from the expansion of the collision term up to second-order. 
The collision term is given by 

^o' = m I iJm^) m<i) iJiw '^''''^•'■' ^'-"'-"'^ i"^!' 

X {s(q')/(p')|l + /(p)|-9(q)/(p)[l + /(p')l} (4.2) 



where ii^(q) ~ {q'^+m'^y^'^, M is the amplitude of the scattering process, S'^{q+p—q'—p') = S^{q+p—q' — p')S{E{q) + 
p — E{c() —p') ensures the energy- momentum conservation and g is the distribution function for electrons. The Pauli 
suppression factors (1 — g) have been dropped since for the epoch of interest the density of electrons rie is low. The 
electrons are kept in thermal equilibrium by the Coulomb interactions with protons and they are non-relativistic, thus 
we can take a Boltzmann distribution around some bulk velocity v 



5(q) = "-e 
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(4.3) 



By using the three dimensional delta function the energy transfer is given by -E(q) — -^(q + p — p') and it turns out 
to be small compared to the typical thermal energies 



E{ci) - i?(q + p - p') ^ ^P P'^ "^ = 0{Tq/m,) . 



(4.4) 



In Eq. (4.4) we have used E{q) = me-|-g^/2me and the fact that since the scattering is almost elastic [p ~ p') (p — p') 
is of order p r^ T, with q much bigger than (p — p'). In general, the electron momentum has two contributions, the 
bulk velocity {q = niev) and the thermal motion (q ^ {raJ'Y''^) and thus the parameter expansion q/nie includes 
the small bulk velocity v and the the ratio (T/me)^/^ which is small because the electrons are non-relativistic. 

The expansion of all the quantities entering the collision term in the energy transfer parameter and the integration 
over the momenta q and q' is described in details in Ref. [31]. Here we just note that it is easy to realize that 
we just need the scattering amplitude up to first order since at zero order g(q') = g(q + p — p') = (7(q) and 
6{E{ci) +p— E{q') — p') = S{p —p') so that all the zero-order quantities multiplying \M\ vanish. To first order 
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where cos^ = n • n' is the scattering angle and aT the Thompson cross section. The resulting collision term up to 
second order is given by [31] 
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where we arrange the different contributions following Rcf. [31]. The first order term reads 

c(i)(p, p') = (1 + cos^^) S{p - p'Kf^'\p') - /(i)(p)) + (/(")(p') - /(")(p))(p - p') • v^fc^ 



(4.6) 



(4.7) 



while the second-order terms have been separated into four parts. There is the so-called anisotropy suppression term 

(4.8) 



4'^(P,P') = ^ (1 + cos^e) 6ip-p'){f^^Hp') /(2)(p)) ; 



a term which depends on the second-order velocity perturbation defined by the expansion of the bulk fiow as v 

v(i)+v(2)/2 



c?Hp, p') = ^(1 + cos2e)(/(°)(p') - /(")(p))(p - p') . v(2) ^^^^^ ; 
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a set of terms coupling the photon perturbation to the velocity 
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and a set of source terms quadratic in the velocity 
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1 + cos^ 



(P ^ P') ■ V d'^S{p - p') 
2 djP 



2cos6'(l-cos6')(n + n') -v 



d5{p — p') 
dp' 



(4.11) 



The last contribution arc the Kompaneets terms describing spectral distortions to the CMB 



„(2) 
-K 



(P.P') ^ (i + ™^.)<E_i^[(/<«.(pO-/'»>(p))r.5!fci2 



/<"'(?') + l"\v) + 2/<»'(p')/(»>(p)) ^fc^ 



2(p — p')coa9{\ — cos^ 



<5b-y)/^"'(p')(l + /^"'(rt) (/("'b') -/^"'(P)) ^^^1^ 



(4.12) 



Let us make a couple of comments about the various contributions to the collision term. First, notice the term 
Cv (p, p') due to second-order perturbations in the velocity of electrons which is absent in Ref. [31]. In standard 
cosmological scenarios (like inflation) vector perturbations are not generated at linear order, so that linear velocities 
are irrotational u'^^' = d'^v^^'. However at second-order vector perturbations are generated after horizon crossing as 
non-linear combinations of primordial scalar modes. Thus we must take into account also a transverse (divergence 
free) component, w'^'* — d'^v^^i -\-Vj, with diVj, — 0. As we will see such vector perturbations will break azimuthal 
symmetry of the collision term with respect to a given mode k, which instead usually holds at linear order. Secondly, 
notice that the number density of electrons appearing in Eq. (4.6) must be expanded as rig = ne(l -|- 6e) and then 



e)c(i)(p,p') 



(4.13) 



gives rise to second-order contributions in addition to the list above, where we split S^ = Se + Se /2 into a first- and 
second-order part. In particular the combination with the term proportional to v in c'^'(p,p') gives rise to the so 
called Vishniac effect as discussed in Ref. [311. 



B. Computation of different contributions to the collision term 



In the integral (4.6) over the momentum p' the first-order term gives the usual collision term 

a/(o) 



C(i)(p)-n,aT 
where one uses the decomposition in Legendre polynomials 



f^'^(p) + o/2''^2(v . n) - /(I) - p-^v ■ n 



dp 



/«(x,p,n) = ^(2^+l)/,(i'(p)P,(cos^), 



(4.14) 



(4.15) 



where t? is the polar angle of n, cosi? = n • v. 

In the following we compute the second-order collision term separately for the different contributions, using the 
notation C(p) = C^^^(p) -I- C^^''(p)/2. We have not reported the details of the calculation of the first-order term 

because for its second-order analog, c^'^(p,p') + cj, (p,p'), the procedure is the same. The important difference is 
that the second-order velocity term includes a vector part, and this leads to a generic angular decomposition of the 
distribution function (for simplicity drop the time dependence) 



/^'Hx,P,n) 



EE 



/«(x,p)(-.)^ / 



4tt 
2i+l 



Ylmin). 



(4.16) 



such that 



if = (-0- 



2^+1 
47r 



df7/«y/„(n) 



(4.17) 



Such a decomposition holds also in Fourier space [34]. The notation at this stage is a bit confusing, so let us restate 
it: superscripts denote the order of the perturbation; the subscripts refer to the moments of the distribution. Indeed 
at first order one can drop the dependence on m setting to = using the fact that the distribution function does not 

depend on the azimuthal angle 0. In this case the relation with /j is 



/l™=H)"'(2^+l)'5™o/f . 



(4.18) 
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a) c'i\p,p') 

The integral over p' yields 



^(.) (p) . ^ I dp'p'^c^l^ (p, p') = ^ / dpWp - p') / ^ [1 + (n • n')^] [/ (^) (P') - /(^) (P)] . 

(4.19) 

To do the angular integral we write the angular dependence on the scattering angle cos0 — n • n' in terms of the 
Legendre polynomials 



[1 + (n • n' 



l\2] 



l + lp2(n.n') 



1 + ^E Wn)r;„>')^ 



where in the last step we used the addition theorem of spherical harmonics 



Pf 



47r 



2e- 



- J2 YUn)Y;^{n') 



(4.20) 



(4.21) 



m=~2 



Using the decomposition (4.17) and the orthonormality of the spherical harmonics we find 



.(2) 



Ca (P) = "-e^T 



f^^iP)-f''\p)-'^ E i^f^'2ip)Y2M 



2 '^ 53/2 

rn— — 2 



(4.22) 



It is easy to recover the result for the corresponding first-order contribution in Eq. (4.14) by using Eq. (4.18). 

b) c('^(p,p') 

Let us fix for simplicity our coordinates such that the polar angle of n' is defined by fi' — v"-^^ • n' with cp' 

(2) 
the corresponding azimuthal angle. The contribution of Cy (p,p') to the collision term is then 

o /• fi^f f\ /*! J / f27v 111 



pV^"Hp') f'''ip)f-^^^ j ^ ^(PA^ - pV) ^ ' ^[1 + (p • P')^] . (4.23) 



We can use Eq. (4.20) which in our coordinate system reads 



[1 + (n • n') 



/n2i 



l + ip2(n-n') 



-J m 



(2-m)! 



2 ^ f2 + m)! 

m— — 2 ^ ' 



P2"(n-v(2^)P2"(n'-v 



-(2)n im(0'- 



so that 



y ^P2(n • n') = P2(n • v(^))P2(n' • v^^)) = P2(/i)P2(A^') • 
By using the orthonormality of the Legendre polynomials and integrating by parts over p' we find 



cfKp) 



Ofio) 

-UeaxP^ — v^-' -n. 
op 



.(2) . . 



(4.24) 



(4.25) 



(4.26) 



As it is clear by the presence of the scalar product v'^) . p ^i^^ final result is independent of the coordinates choosen. 

c) 4'i(p.p') 

Let us consider the contribution from the first term 

4'i(/)(P.P') - (1 + cos^^) (/(^)(P') - f^'Hp)) (P - P') • v^fc^ , (4.27) 
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where the velocity has to be considered at first order. In the integral (4.6) it brings 



:C' 



(2) 



SngCTTW 



2"'A«(/) 4p 






p'^.')fj ^{l+cosH 



(4.28) 



The procedure to do the integral is the same as above. We use the same relations as in Eqs. (4.24) and (4.25) where 
now the angles are those taken with respect to the first-order velocity. This eliminates the integral over (p' , and 
integrating by parts over p' yields 



9^At,(7)(P) - 



3ne(TTV f dji' '" 
4p 7-1 T 



\ + \P2{p)P2{l^') 



Pi^, - 2M')(.f(^)(p,M') - ./(^Hp,/^)) +P^(M - m')-^^^'^^'''^'^ 



dp 



We now use the decomposition (4.15) and the orthonormality of the Legendre polynomials to find 



(4.29) 



/^m7(^'(p,m')/^2(//) = Y.j'^4^'P^(^')P^(^^')f'"'^P)-Y.j''-Y 



^Pi(^') + ^^3(m') 



Pi{^^').ff\p) 



lA'\p) + lfi'\p). 



where we have used n'P2{n')Pi{n') = |Pi(a*') + ^Psip'), with Pi(m') = m'- Thus from Eq. (4.29) we get 



.(2) 



o^AvmiP) = neCrT<V-n 



'f^'Hp) - f!>'\p)-p'-^ - ip.(v . n) (f^'\p)+p'ji^^ 



dp 



+ V 



(1)^ 



(l).^^ , Jfi'iP) . 1 



2,/r^(p)+p 



^^ +-P2(v-n)l2/}^^(p)+p 



9p 



a, +3/-^(,) + -p- ^^ 



In c'-^-'(p, p') there is a second term 



(4.30) 



.(4.31) 



„(2) 



-A.(/7) = -2cos0(l - COS0) [.r\p') - r'\p)) 5{p - p'){n + n') • v , 
whose contribution to the collision term is 



2^Av{II) 



(2) 



(P) 



Suearv 



.(2) 



(4.32) 



j dp'p'5(p-p')j ^{fW{p')-fW{p)){^ + ^')l^^cos0{l-cos9). (4.33) 



This integration proceeds through the same steps as for C/^Lj)ip)- In particular by noting that cos^(l — cos0) = 
-1/3 + Pi(cos0) - 2P3(cos0)/3, Eqs. (4.24) and (4.25) allows to compute 



-cos0(l — cos6 



27r 



and using the decomposition (4.15) we arrive at 



^ + P,{p)P,{f,')-^P,{p)P,{t,') 



1 



.(2) 



(l)r 



2^A.(//)(P) - -^eCTT |v . nf^'>{p){l - P2(v • n)) + 
Thus summing Eq. (4.29) and (4.35) we obtain 
1 



ip2(v-n) (3/«(p)-3/P(p) 



(4.34) 



(4.35) 



7;CaI(p) = ne(7T<v-n 



p(i) 



f^'Hp)-.4'Hp)-p'-^^^^-^''' 



Op -f^"(p) + lp^(^-)[f^''ip)-p^^ 



Vi%)-P'-^ + ^-2(. . n) (-/-(.) ..^ . 6/i^)(.) . yj^ 



.(4.36) 



As far as the remaining terms, these have already been computed in Ref. [31] (see also Ref. [30]) and here we just 
report them 
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d) 



c™(P.P') 



The term proportional to the velocity squared yield a contribution to the collision term 



77C™ (P) ^nearl (v • n 



9/(0) 11 2^^^ 



9/(0) 
dp 



,Oy(o) 



20^ 9p2 



(4.37) 



e) c(^^(p,p') 

The terms responsible for the spectral distortions give 



op 



Finally we write also that part of the collision term coming from Eq. (4.13) 



e^ cW(p,p') - <5WcW(p) = ne^T^W 



eWi 



f(i 



/rb)+o/2 ^2(v-n)-/ 



(1) 



P- 



9/W 
dp 



(4.38) 



(4.39) 



C. Final expression for the collision term 

Summing all the terms we find the final expression for the collision term (4.6) up to second order 

C(p) = C«(p) + ic(2)(p) 



with 



and 



C^'\p)^n,aT 



fW, 



f(i) 



/rW+0/2 ^2(v-n)-/ 



_ f(i) _ 



a/(o)_ 
dp 



-V • n 



ic(^)(p) = n.ar\lCip)-\ E ^ /^^ W ^2™(n) - i/(^)(p) + ^ 
— —p^ — v*- -' • n + V • n 



fWi 



f(i) 



(4.40) 



(4.41) 



/r(p)+o/2 ^2(v-n)-/ 



(1) 



dp 



dp 



f^'Hp)-f^'\p)-P^^^-^''^ 



w- " -^"'^^^ + ^"^^^^ • "^ (-^^'''^^^ -^^€^ 



2/} nrt +P gp + 5^2(v • n) I -/} ^(p) +P g^ + 6/;^ ^(p) + -p "g^ 



+ (v-n)2 



a/(") n 2 9'/^°^ 

' dp ^ 20^ 9p2 



9p ' 20^^ 9p 



7^P 



,a2j(0) 






pM7^e^+.f(°)(l + ./(°)) 



(4.42) 



Notice that there is an internal hierarchy, with terms which do not depend on the baryon velocity v, terms 
proportional to v • n and then to (v • n)^, v and v^ (apart from the Kompaneets terms). In particular notice the term 

proportional to 5e v • n is the one corresponding to the Vishniac effect. We point out that we have kept all the terms 
up to second-order in the collision term. In Refs. [30, 31] many terms coming from c)^^ have been dropped mainly 
because these terms are proportional to the photon distribution function f^^> which on very small scales (those of 
interest for reionization) is suppressed by the diffusion damping. Here we want to be completely general and we have 
to keep them. 
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V. THE BRIGHTNESS EQUATION 



First-order 



The Boltzmann equation for photons is obtained by combining Eq. (3.24) with Eqs. (4.41)-(4.42). At first-order the 
left-hand side reads 



4f 
drj 



dfw a/(") a$(i) dx' a/(o) 9^(1) 



dr] 



P- 



dp dx^ drj 



P- 



dp drj 



(5.1) 



At first-order it is useful to characterize the perturbations to the Bose-Einstein distribution function (3.22) in terms 
of a perturbation to the temperature as 



f{x\p,n\i]) = 2 



exp 



P 



Thus it turns out that 



/ 



(1) 



r(77)(i + e(i)) 
a/(o) 



- 1 



dp 



.6(1) 



(5.2) 



(5.3) 



where we have used that 9//99|e=o = —pdf^^'' /dp. In terms of this variable O^^^ the linear collision term (4.41) will 
now become proportional to —pdj^^' /dp which contains the only explicit dependence on p, and the same happens 
for the left-hand side, Eq. (5.1). This is telling us that at first order O'^) does not depend on p but only on x*, n', 77, 
Q(i) = ©(1) (x% n*, r). This is well known and the physical reason is that at linear order there is no energy transfer in 
Compton collisions between photons and electrons. Therefore the Boltzmann equation for Q^^> reads 



aew 

drj 



BQW a$(i) 



dx'^ 



dx^ 



drj 



n^UTCL 



e 



(1) 



ie«P2(v-n) 



Qd) 



(5.4) 



where we used that /^ = —pdf^^'' /dpQ\ according the decomposition of Eq. (4.15), and we have taken the zero- 
order expressions for dx'' /drj, dropping the contribution from dn^ /drj in Eq. (3.8) since it is already first-order. 
Notice that since &^^' is independent on p it is equivalent to consider the quantity 



A'^^\x\n\T) ^ 
A(i) = 



Jdpp'^f^) 

J dpp3f(0) ' 

49(1). 



(5.5) 
(5.6) 



The physical meaning of A(i) is that of a fractional energy perturbation (in a given direction). From Eq. (3.24) 
another way to write an equation for A^^' - the so called brightness equation - is 



d 

drj 


A(i)+4$(i)' 


_ 4|- ($(1) + vl/d)) 


= neOTl 



A 



(1) 



A^^'P2(v-n)-A(i)+4vn 



(5.7) 



B. Second-order 



The previous results show that at linear order the photon distribution function has a Planck spectrum with the 
temperature that at any point depends on the photon direction. At second-order one could characterize the perturbed 
photon distribution function in a similar way as in Eq. (5.2) 



f{x\p,n\Tj) = 2 



exp 



P 



T{rj) e^- 



- 1 



n -1 



(5.8) 



where by expanding O = 0(i) + 6(^)/2 + ... as usual one recovers the first-order expression. For example, in terms of 
0, the perturbation of f^^> is given by Eq. (5.3), while at second-order 



/(2) 



2 dp 2 



,^2/(0) a/(0) 



dp"^ 



dp 



0(1) 



(5.9) 
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However, as discussed in details in Refs. [30, 31], now the second-order perturbation Q^"^^ will not be momentum 
independent because the collision term in the equation for O'^' does depend explicitly on p (defining the combination 
— (p9/(°V^P)~^/^^^ does not lead to a second-order momentum independent equation as above). Such dependence is 
evident, for example, in the terms of C'^'(p), Eq. (4.42), proportional to v or v"^, and in the Kompaneets terms. The 
physical reason is that at the non- linear level photons and electrons do exchange energy during Compton collisions. 
As a consequence spectral distorsions are generated. For example, in the isotropic limit, only the Kompaneets terms 
survive giving rise to the Sunyaev-Zeldovich distorsions. As discussed in Ref. [30], the Sunyaev-Zeldovich distorsions 
can also be obtained with the correct coefficients by replacing the average over the direction electron (w^) with the mean 
squared thermal velocity (v^;j) = 3Te/me in Eq. (4.42). This is due simply to the fact that the distinction between 
thermal and bulk velocity of the electrons is just for convenience. This fact also shows that spectral distorsions due 
to the bulk flow (kinetic Sunyaev-Zeldovich) has the same form as the thermal effect. Thus spectral distorsions can 
be in general described by a global Compton y-parameter (see Ref. [30] for a full discussion of spectral distorsions) . 
However in the following we will not be interested in the frequency dependence but only in the anisotropics of the 
radiation distribution. Therefore we can integrate over the momentum p and define [30, 31] 



A(2)(a 



3 f (2) 



J dppSf^o) ' 



(5.10) 



as in Eq. (5.5). 

Integration over p of Eqs. (3.24)-(4.42) is straightforward using the following relations 



dpp p 



3 5/(") 



dpp p 



dp 

3 2aV(°) 



-AN . 



20 7V, 



dp'^ 
dpp^f^^^ = 7VA(i) 



dpp p 



dfW 



dp 



= -47VA(i) 



r,3 f (0) 



N = / dpp'f 



(5.11) 



where N is the normalization factor in Eq. (5.10) (it is just proportional the background energy density of photons 
p-y). At first-order one recovers Eq. (5.7). At second-order we find 



1 d r 
2dr] 



A(2) +4$(2) 



d 

drj 



U(i)+4$« 



4A(i)f*(i)'-$(i), 



2I- U(^) + $(2)) + 4^n^ + 2^nV^- 
drj \ J or] or) 



T 
"2" 



^IS - A(^) - ^ E tS ^-n y2m{n) + 25(^) ( A« + iA(^)p,(v . n) - A^^^ + 4v • n ) + 4v(^) • „ 



2 ^53/2 

m— — 2 



-^2(v • n) 



A(i)-h3A 



(1) a(i) 



2 (l--P2(v-n) 



- vA[^'^ (4 + 2P2(v • n)) + 14(v • n)^ - 2w^ 



(5.12) 



where we have expanded the angular dependence of A as in Eq. (4.16) 



A«(x,„) = E E <lw( 



4tt 



2i+l 



y«m(n). 



(5.13) 



with 



A 






HY 



2i+l 
47r 



df7A(*)y;,„(n) , 



(5.14) 



where we recall that the superscript stands by the order of the perturbation. At first order one can drop the dependence 
on m setting m == so that A^^ = {—i)^^{2£ + l)Smo ^e ■ In Eq- (5-12) we have introduced the differential optical 
depth 



-UelTa . 



(5.15) 
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It is understood that on the left-hand side of Eq. (5.12) one has to pick up for the total time derivatives only those 
terms which contribute to second-order. Thus we have to take 



1 d 
2dr) 



A(2) +4$(2) 



d 

drj 



A(i)+4$« 



(2) _ If d_ ^_d_ 
2 \dri dx"^ 



($W + ^W)n*n^' - ($■* -F *'*) 



A(2) + 4$(2)) + „*($(!) + $(i))a^(A(i) + 4$(i)) 



9n* 



(5.16) 



where we used Eqs. (3.10) and (3.21). Notice that we can write aA(i)/9n* = {dA^'^"> /dx'){dx' /dn') = {d/S.^'^"> /^x^){rJ- 
^q^)^ from the integration in time of Eq. (3.10) at zero-order when n* is constant in time. 

C. Hierarchy equations for multipole moments 

Let us now move to Fourier space. In the following for a given k mode we will choose the coordinate system such that 
63 = k and the polar angle of the photon momentum is -d, with /i = cosz? = k • n. Then Eq. (5.12) can be written as 

A(2)' 4- z/c^A(2) _ r'A(2) ^ 5'(k, n, 77) , (5.17) 

where S'(k, n, 77) can be easily read off Eq. (5.12). We now expand the temperature anisotropy in the multipole 

moments A)^ in order to obta 
Eq. (5.14) to Eq. (5.17) we find 



moments A)^ in order to obtain a system of coupled differential equations. By applying the angular integral of 



A£!'(k,ry) = A: 



«^m . (2) «£+l,m . (2) 

2^-1 ^-1'™ ^2^+3 ^+1'™ 



T'/\f}^ + Sim , (5.18) 



where the expansion coefficients of the source term are given by 



Stm - (4*^''' - r'A^'o^) 5n6mo + 4fc$(2)j^^^^p _ 4cj±i(5,i - 8r'«^),5,i - ^A^^ fe - 2x±2 fe 



2r 



, f d^ki 



vl^H^iK\^2)^2 ■ ki + 5(i)(ki)A^^^(k2) - ^-vik,)A[l\k2 



(27r)3 



(27r)3 
, f d^ki 



SeiSmo — 2 



(vl/(i)v$(i)),„ + 8t'((5(i)v)„, + +6t'(A(^'v),„ - 2t'(A(^V), 



5n 



(2^) 



S(^\k,)Ai'\k2) - ^-v{k,)A[l\k,) 



Sn <5mO + 



d^h 

'(2nf 



8*(i)(ki) + 2T'(5(i)(ki) 



iv - V^)i'^^'^ + <^^'^){^i) ki • ksl A(;'(k2)<5„ 



i{-i)-^{-l)-"'(2£ + 1) 5] ^ {2£" + 1) [8A(^?V$(i) - 2($(i) + *(i))VA^,V 



f 1 ^ \ / f 1 £ 
,' V / I m' -m 



m' — — 1 
1 



r'*(-z)-^(-l)-"(2^+l)^ ^ (2r + l) 



£" m' = -l 



.(1) 



2A,Vv + 5(5£"2A^'V 



(1), 



^" 1 ^ \ / ^" 1 e 

V I j I m' -TO 



+ 14t'(-z)-^(-1)-™ 
1 1 £ 



d^ki 

^ J Wyi 

1 1 i 



E 



(1)^ 



„a,„,,1^4. 



yi;,.(i<i) a'i;."(k)-*'ir,:„.,(ki) 



/ I TO,' ?7i" —ni 
+ 2(r7-ry,)(-^)-^(-l)- 



2^+1 



E E 



47r ^^ ^^ 7 (27r)3 V 2L + 1 V 3 

L m' .m" — — 1 



''''^ ^^^^^VA(^)(kO($«+*(^))(k2) 



X fciF, 



i*„,(ki) (fcri:„„(k) - fciYi:„„(ki)) y'rfr!yi„,(n)yi„»(n)rio(n)r,-™(n) , 



(5.19) 



where 



k2 = k ki , 



(5.20) 



16 



and ^2 — |k2|. In Eq. (5.19) it is understood that \m\ < L 

Let us explain the notations we have adopted in writing Eq. (5.19). The baryon velocity at linear order is irrotational, 
meaning that it is the gradient of a potential, and thus in Fourier space it is parallel to k, and following the conventions 
of Rcf . [37] , we write 



v(i)(k) = -z4^^(k)k. 



(5.21) 



For the second-order velocity perturbation, it will contain a transverse (divergence-free) part whose components are 
orthogonal to k = 63, and we can write 



v(^)(k) = -..f(k)e3+ E-L^'^^> 



m=±l 



(5.22) 



where e^ form an orhtonormal basis with k. The second-order perturbation uoi is decomposed in a similar way, and 
LO±i are the corresponding components (in this case in the Poisson gauge there is no scalar component). Similarly for 
the tensor perturbation Xij we have indicated its amplitudes as x±2 in the decomposition [36] 



Xi] 




(ei±ie2)i(ei ± ie2) 



(5.23) 



We have taken into account that in the gravitational part of the Boltzniann equation and in the coUsion term there 
are some terms, like 5e v, which still can be decomposed in the scalar and transverse parts in Fourier space as in 
Eq. (5.22). For a generic quantity /(x)v we have indicated the corresponding scalar and vortical components with 
(/v)m and their explicit expression is easily found by projecting the Fourier modes of /(x)v along the k = 63 and 
(e2 =F *Gi) directions 



(/v)™(k) 



^««(ki)/(k2)yr„(ki)\/^. 



Similarly for a term like /(x)V(7(x) we used the notation 



(/V5)™(k) = - 



-fcig(ki)/(k2)Yi:„(ki)\/— . 



(27r)3 '"^ ^'■"- '■' ""^ ^'V 3 
Finally, the first term on the right-hand side of Eq. (5.18) has been obtained by using the relation 



^k.nA(2)(k) = EAf2(k) 



2^+1 



HemGi-i^rn ~ K^+l.mGf+i^m 



^E 



l^em » (2) _ l^im » (2) 

2^-1 ^-i'" 2^ + 3 ^+1'™ 



where Gim — (— 2)^-\/47r/(2i' + l)Yfm(n) is the angular mode for the decomposition (5.13) and 



Kir, 



V^ 



(5.24) 
(5.25) 

(5.26) 
(5.27) 



This relation has been discussed in Refs. [34, 36] and corresponds to the term n^SA^^^/Qa;* in Eq. (5.12). 

As expected at second-order we recover some intrinsic effects which are characteristic of the linear regime. In 
Eq. (5.18) the relation (5.26) represents the free streeming effect: when the radiation free streems the inhomogeneities 
of the photon distribution are seen by the obsever as angular anisotropics. At first order it is responsible for the 
hierarchy of Boltzmann equations coupling the different i modes, and it represents a projection effect of fiuctuations 
on a scale k onto the angular scale £ ~ kr]. The term t'A)^ causes an exponential suppression of anisotropics in the 
absence of the source term Sim- The first line of the source term (5.19) just reproduces the expression of the first 
order case. Of course the dynamics of the second-order metric and baryon-velocity perturbations which appear will 
be different and governed by the second-order Einstein equations and continuity equations. The remaining terms in 
the source are second-order effects generated as non-linear combinations of the primordial (first-order) perturbations. 
We have ordered them according to the increasing number of £ modes they contribute to. Notice in particular that 
they involve the first-order anisotropics A^ and as a consequence such terms contribute to generate the hierarchy 
of equations (apart from the free streeming effect). The source term contains additional scattering processes and 
gravitational effects. On large scales (above the horizon at recombination) we can say that the main effects are due 
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to gravity, and they include the Sachs- Wolfe and the (late and early) Sachs- Wolfe effect due to the redshift photons 
suffer when travelling through the second-order gravitational potentials. These, toghether with the contribution due 
to the second-order tensor modes, have been already studied in details in Ref. [19]. Another important gravitational 
effect is that of lensing of photons as they travel from the last scattering surface to us. A contribution of this type 
is given by the last term of Eq. (5.19). On small scales scattering effects around the epoch of recombination are 
important and we will study them in details in a companion paper [27] . 

D. Integral solution to the second-order Boltzmann equation 

As in linear theory, one can derive an integral solution to the Boltzmann equation (5.12) in terms of the source term 
S. Following the standard procedure (see e.g. Ref. [28, 35]) for linear perturbations, we write the left-hand side as 
A(2)' + ikfiM"^^ - r'A(2) = e-*'=''+^d[A(2)e*'=^''-^]/d77 in order to derive the integral solution 



rva 
A(2)(k,n,ryo)= / d7?5(k,n,r?)e*'=^^(''-'"')-^, (5.28) 

(2) 

:ime. The expression of the photon moments A^^ can be obtained as usual fri 
Eq. (5.14). In the previous section we have already found the coefficients for the decomposition of source term 5* 



(2) 

where r/o stands by the present time. The expression of the photon moments A^^ can be obtained as usual from 



^(k,n,ry) ^J2J2 Sem{Kv){-^f J ^l^Y^n) . (5.29) 

In Eq. (5.28) there is an additional angular dependence in the exponential. It is easy to take it into account by 
recalling that 

g.k.x ^ ^(j)^(2£ + l)j,{kx)Pe{k ■ x) . (5.30) 

e 
Thus the angular integral (5.14) is computed by using the decomposition of the source term (5.29) and Eq. (5.30) 

Af2(k,%) - (-l)-"(-z)-^(2^+l) / rfr^e--^")^ ^ {-zY^S,,^,J2'''{2e^ + l)j,,{k{v - Vo)K2ii + 1) 

^ [o j (^ m'2 ~m j ' (^-^^^ 

where the Wigner 3 — j symbols appear because of the Gaunt integrals 

{2h + 1) (2^2 + 1) (2^3 + 1) ( h i2 h\ ( il i2 is 



4tt V / y TTii 1712 "m-s 

Since the second of the Wigner 3-j symbols in Eq. (5.31) is nonzero only if to = to.2, our solution can be rewritten to 
recover the corresponding expression found for linear anisotropics in Refs. [34, 36] 

^^^=/ rfr^e-C')^ Y. S,,^J--\ki,jo-v)], (5.32) 

where j) ^'"^ [k{rio — rj)] are the so called radial functions. Of course the main information at second-order is included 
in the source term containing different effects due to the non-linearity of the perturbations. In the total angular 
momentum method of Refs. [34, 36] Eq. (5.32) is interpreted just as the intergration over the radial coordinate 
(x — % ~ '7) of the projected source term. Another important comment is that, as in linear theory, the integral 
solution (5.31) is in fact just a formal solution, since the source term S contains itself the second-order photon 
moments up to / = 2 (see Eq. (5.19)). This means that one has anyway to resort to the hierarchy equations for 
photons, Eq. (5.18), to solve for these moments. Neverthless, as in linear theory [35], one expects to need just a 
few moments beyond ^ := 2 in the hierarchy equations, and once the moments entering in the source function are 
computed the higher moments are obtained from the integral solution. Thus the integral solution should in fact be 
more advantageous than solving the system of coupled equations (5.18). 
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VI. THE BOLTZMANN EQUATION FOR BARYONS AND CDM 



In this section we will derive the Boltzniann equation for massive particles, which is the case of interest for baryons 
and dark matter. These equations are necessary to find the time evolution of number densities and velocities of the 
baryon fluid which appear in the brightness equation, thus allowing to close the system of equations. Let us start from 
the baryon component. Electrons are tightly coupled to protons via Coulomb interactions. This forces the relative 
energy density contrasts and the velocities to a common value, S,^ = Sp = Sf, and v^ = Vp = v, so that we can identify 
electrons and protons collectively as "baryonic" matter. 

To derive the Boltzmann equation for baryons let us first focus on the collisionless equation and compute therefore 
dg/drj, where g is the distribution function for a massive species with mass m. One of the differences with respect 
to photons is just that baryons are non-relativistic for the epochs of interest. Thus the first step is to generalize the 
formulae in Section II up to Eq. (3.21) to the case of a massive particle. In this case one enforces the constraint 
Q^ = g^ivQ^Q^ = —m? and it also useful to use the particle energy 



E ~ \J cp' + rn? . 



(6.1) 



where q is defined as in Eq. (3.3). Moreover in this case it is very convenient to take the distribution function as a 
function of the variables (f — qv} , the position x* and time r], without using the explicit splitting into the magnitude 
of the momentum q (or the energy E) and its direction n*. Thus the total time derivative of the distribution functions 
reads 



(6.2) 



dg dg dg dx^ dg dq^ 
drj drj dx"^ drj dq^ drj 

We will not give the details of the calculation since we just need to replicate the same computation we did for the 
photons. For the four- momentum of the particle notice that Q* has the same form as Eq. (3.6), while for Q^ we find 



,-* 



■Ell + u;,^ 



(6.3) 



In the following we give the expressions for dx^ /drj and dq^ /drj. 
a) As in Eq. (3.10) dx^ /d-q — Q^ /Q^ and it turns out to be 



— — — — n e ^ 
d-q E 



^'"* f) ( ^ ^ 'KXtmn'^n''' 



(6.4) 



b) For dq^ /drj we need the expression of Q* which is the same as that of Eq. (3.6) 



a \ Z 



(6.5) 



The spatial component of the geodesic equation up to second-order reads 

f = -2(H-*')fl-Jx.™nV")^n^e* + 2|^4nVe*-*-|i^e 
drj \ 2 J a ox^ aE ox^ a 



$+2<I' 



E 



a 



{x\ + ^'fe - O- + n^'5jk - [x'jM + x\j + x.,k) 



q^q^ 
Ea 



a* q^ 

dx^ aE 



*+2<I' 



{u'' + Hlo') 



E 



(6.6) 



Proceeding as in the massless case we now take the total time derivative of Eq. (6.5) and using Eq. (6.6) we find 

~d^ '■- ^ '-^ ' ^'- E " ^ "" ^''S' 



= -{n- ■^')q' + vl/,fci|-e*+* - $'^Se*+* - *.5e*+* - E{lo'' + 'H^l - {x"k + ^'k - OE 



Hij'S.k - {x\.k + X'fcj + X,fe) 



q^c^ 



(6.7) 



We can now write the total time derivative of the distribution function as 



drj drj E 



1 — dJiTl^ 



l:Xkmn'nA |^ 
2^ / dx' 



(W -*')?* + *, 






<P'Ee 



$+iii 



*.,;te*+* 



E 



E{u;^' + Hu;^) - {x\ + oj)^ - u;'^)E + (Uuo^d.k - (x^-,fe + x\,, + X.O) 



q^q'^ 



55 
dq^ 



(6.8) 
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This equation is completely general since we have just solved for the kinematics of massive particles. As far as the 
collision terms are concerned, for the system of electrons and protons we consider the Coulomb scattering processes 
between the electrons and protons and the Compton scatterings between photons and electrons 

-^(X,q,?7) = {Cep)QQ'q' + {Ce^)pp'q' (6.9) 

-^(X, Q,?7) == {Cep)qq'Q' , (6.10) 

where we have adopted the same formalism of Ref. [28] with p and p' the initial and final momenta of the photons, q 
and q' the corresponding quantities for the electrons and for protons Q and Q'. The integral over different momenta 
is indicated by 

/ \ = f ^L. [ fjL [ ^^ ((^^^^ 

{■■■)pp','-j (2^)3 y (2^)3 y (2^)3 •••' ^^■'') 

and thus one can read c^'y as the unintegrated part of Eq. (4.2), and similarly for Cep (with the appropriate amplitude 
|Mp). In Eq. (6.9) Compton scatterings between protons and photons can be safely neglected because the amplitude 
of this process has a much smaller amplitude than Compton scatterings with electrons being weighted by the inverse 
squared mass of the particles. 

At this point for the photons we considered the perturbations around the zero-order Bose-Einstein distribution 
function (which arc the unknown quantities). For the electrons (and protons) we can take the thermal distribution 
described by Eq. (4.3). Moreover we will take the moments of Eqs. (6.9)-(6.10) in order to find the energy-momentum 
continuity equations. 

A. Energy continuity equation 

We now integrate Eq. (6.8) over d^q/{2'!r)^. Let us recall that in terms of the distribution function the number density 
Ue and the bulk velocity v are given by 

^g, (6.12) 



and 



(27r) 



1 f (Pc 



qn 



Tie J (27r)3 " E 



(6.13) 



where one can set E ~ rUe since we are considering non-rclativistic particles. We will also make use of the following 
relations when integrating over the solid angle (Kl 

drtn' ^ I dn n'n^n^ = , / — n^n> ^ -6'^ . (6.14) 

J J 47r 3 

Finally notice that dE/dq = q/E and dg/dq — {q/E)dg/dE. 

Thus the first two integrals just brings n^ and (neW*),i. Notice that all the terms proportional to the second-order 
vector and tensor perturbations of the metric give a vanishing contribution at second-order since in this case we can 
take the zero-order distribution functions which depends only on j] and E, integrate over the direction and use the 
fact that S^^Xij = 0- The trick to solve the remaining integrals is an integration by parts over <;*. We have an integral 
like (the one multiplying (^' — H)) 

d^q ,■ dg f d^q 

(2^'^a^ = -^y (2^^ = -^"- (^-^^^ 

after an integration by parts over q^. The remaining integrals can be solved still by integrating by parts over q*. The 
integral proportional to <I>^* in Eq. (6.8) gives 



^ *' / (2^^# - ^ ^'"- (6-16) 
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where we have used that dE/dq^ = <t /E. For the integral 

the integration by parts brings two pieces, one from the derivation of q^q^ and one from the derivation of the energy 
E 

- 4e*-** . / 4^A + e---^., I 4^yii - -4e*+**,..^ + e^--^, f ^^g^-^ . (6.18) 



(27r)3-^£; '"y {271^ E E ' '" ■] {2t:Y'' E"^ E 

The last integral in Eq. (6.18) can indeed be neglected. To check this one makes use of the explicit expression (4.3) 
for the distribution function g to derive 

dg qt me 

g— -77r«»5, (6.19) 



and 



Thus it is easy to compute 



d^q 
(27r)3 



dv^ " Te Te 



gq'q^ == S'^fienieTe + nemlv'v^ . (6.20) 



^$+»I>^,fc f d q 2„fc _ ^*+>I',T; ,,2 -^"^ _L q^*+*iT( ,, ^ '^^ ^ ^*+*> 



6 — T / TTT'^m q ^ -e ^ *,fcU h 3e ^ *,fci'fe'^e Ke ^ *,fcWfcW , (6.21) 

m^ J (2tt)-^ me me 

which is negligible taking into account that Te/me is of the order of a thermal velocity squared. 

With these results we are able to compute the left-hand side of the Boltzmann equation (6.9) integrated over 
d^q/{2iTY . The same operation must be done for the collision terms on the right hand side. For example for the first 
of the equations in (6.9) this brings to the integrals {cep)QQ'qq' + {ce-^) pp' qq' ■ However looking at Eq. (4.2) one realizes 
that {ce'y)pp'qq' vanishcs because the integrand is antisymmetric under the change q <-^ q' and p <-^ p'. In fact this 
is simply a consequence of the fact that the electron number is conserved for this process. The same argument holds 
for the other term {cep) qq' qq' ■ Therefore the right-hand side of Eq. (6.9) integrated over d'^q/ {2Tr)^ vanishes and we 
can give the evolution equation for rig. Collecting the results from Eq. (6.15) to (6.21) we find 

dne ^ ^$+3, 9(wVie) ^ 3.^ _ ^,. _ 2g*+*^ ^^k ^ g*+*^ kv'' = 0. (6.22) 

or] ox^ 

Similarly, for CDM particles, we find 

— ^ \-e — h 3(/t - W )ncDM - 2e y^^feWcDM + ^ '±',fcWcDM = L) • (6.23) 

B. Momentum continuity equation 

Let us now multiply Eq. (6.8) by [q'' / E) / {2Tr)^ and integrate over d^q. In this way we will find the continuity equation 
for the momentum of baryons. The first term just gives (rief *)'. The second integral is of the type 

where we have used Eq. (6.20) and E = me- The third term proportional to (7i — ^') is 

d^q kdg q' . ^ f d^q q^ q' 

-4ne+ / TTT-TT.?!^!^, (6.25) 



{2tt)^^ dqkE " J {2nf^E^ E' 

where we have integrated by parts over q*. Notice that the last term in Eq. (6.25) is negligible being the same integral 
we discussewd above in Eq. (6.21). By the same arguments that lead to neglect the term of Eq. (6.21) it is easy to 
check that all the remaining integrals proportional to the gravitational potentials are negligible except for 
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The integrals proportional to the second-order vector and tensor perturbations vanish as vector and tensor pertur- 
bations are traceless and divergence-free. The only one which survives is the term proportional to cj" -I- Tiw* in 
Eq. (6.8). 

Therefore for the integral over S'qq^ /E of the left-hand side of the Boltzmann equation (6.8) for a massive particle 
with mass rrie {'nip) and distribution function (4.3) we find 



4{H - *')"e«' + *^*e*+*ne + e*+* n^^- + e*+*---(neW^w*) + ^n^ + Hlu'ii, 



{2Tr)^ E drj drj \ me J dxi drj 

(6.27) 

Now, in order to derive the momentum conservation equation for baryons, we take the first moment of both 
Eq. (6.9) and (6.10) multiplying them by q and Q respectively and integrating over the momenta. Since previously 
we integrated the left-hand side of these equations over d^qq^ / E, we just need to multiply the previous integrals by 
TOg for the electrons and for irip for the protons. Therefore if we sum the first moment of Eqs. (6.9) and (6.10) the 
dominant contribution on the left-hand side will be that of the protons 



(fiQ J dgp 



(27r)3 ^ dr/ 



^ {Cep{q" + Q'))QQ'qq' + {Ce^q')pp' qq' ■ (6.28) 



Notice that the integral of the Coulomb collision term Cep{(f + Q*) over all momenta vanishes simply because of 
momentum conservation (due to the Dirac function 5'^{q + Q — q' — Q')). As far as the Compton scattering is 
concerned we have that, following Ref. [28], 

\Ce^q ) pp' qq' — {pe^P /pp'qq' 7 yO.ZiJj 

still because of the total momentum conservation. Therefore what we can compute now is the integral over all momenta 
of Ce-yP*. Notice however that this is equivalent just to multiply the Compton collision term C(p) of Eq. (4.2) by p^ 
and integrate over d^p/{2TT^) 



d^p 

(2^ 



{CeyPlpp'qq' = I JTA^fCip) . (6.30) 



where C(p) has been already computed in Eqs. (4.41) and (4.42). 

We will do the integral (6.30) in the following. First let us introduce the definition of the velocity of photons in 
terms of the distribution function 



d^p 

(2^ 



(p^+p>;= / j^fp\ (6.31) 



where p-y — p-y/S is the photon pressure and p^ the energy density. At first-oder we get 

iyW^^I^AWn\ (6.32) 

where A is the photon distribution anisotropics defined in Eq. (5.10). At second-order we instead find 

4 4^^1 frff7 ^,^ ,_4 

3227 47r 3^^ ^ ' 

Therefore the terms in Eqs. (4.41) and (4.42) proportional to /^^■'(p) and /'•^•'(p) will give rise to terms containing 
the velocity of the photons. On the other hand the terms proportional to /g (p) and /gg (p), once integrated, vanish 
because of the integral over the momentum direction n*, j dVtn'^ = 0. Also the integrals involving P2(v • n) = 
[3(v • n)2 - l]/2 in the first fine of Eq. (4.41) and (4.42) vanish since 

/ dl7P2(v ■n)n' ^ v'^v^ f dnukUjii' = , (6.34) 

where we are using the relations (6.14). Similarly all the terms proportional to v, (v • n)^ and v'^ do not give any 
contribution to Eq. (6.30) and, in the second-order collision term, one can check that /dJlF2(n)n* = 0. Then there 
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are terms proportional to (v • n)/*^''^(p), (v • n)pdf^'^^ /dp and (v • lajpOf^ /dp for which wc can use the rules (5.11) 
when integrating over p while the integration over the momentum direction is 



Finally from the second line of Eq. (4.42) we get three integrals. One is 

where p^ is the background energy density of the photons. The second comes from 
^ [ d^P ,, ,^ ,, . ( Ji). , df^^\p)\ 5 

2 7 72^^^^'"^ '^^■"M'^2 (rt-p 



(6.35) 



(6.36) 



dp 



= 7^A« 



dfl 

— 1 
47r 



ivjVkVi I — ri^nPn n — Vj j - — n^n-' 



dn 

— 1 

4tt 



-p^A««% 



(6.37) 



where we have used the rules (5.11), Eq. (6.14) and J{dn/4Tr)n'n'n''n' = {d'^d^' + d'''d'^ + d''d3'')/15. In fact the 
third integral 



exactly cancels the previous one. Summing the various integrals we find 



(6.38) 



J (2^^^P^P " """"^^^^ 



i(v«-v«) 



driA(2) 4v(2) 4 .1,, (1) 



47r 2 



(6.39) 



Eq. (6.39) can be further simplified. Recalling that S^ — Aq we use Eq. (6.33) and notice that 



y'^AW(v.n)n^ = ««n^« + i«^A« 



(6.40) 



where 11*^ is the quantity defined in Eq. (7.5). 

Thus our final expression for the integrated collision term (6.30) reads 



J (2^^^P^^' " "^''^^'' 



A /,,(2)i ,,(2)A A 

(i)n , 4 / u^ ^ v^ \ ^ l^(i)^„,(i)j „(i) 

3 



^iv^'^' - v^;'i + -A^-^] + -A''(-'''' - -'-r''') + ^'K 



3 \ 2 



(6.41) 



We are now able to give the momentum continuity equation by for baryons by combining mpdgp/dr] from Eq. (6.27) 
with the collision term (6.41) 



\ mp 



dr] 

—ricaTO, Pj 



,*+i' 



d 



duj^ 



■^—{pbv^v') + — -pt + Huj'pb 



dx^ 



drj 



4 /«(2). ^{2)^^ 



3 (««^ - v^^)^) + ir—-^^\+ l^w (.(!)' - .w^) + ««n- 



(6.42) 



where ph is the baryon energy density and, as we previously explained, we took into account that to a good 
approximation the electrons do not contribute to the mass of baryons. In the following we will expand explicitly at 
first and second-order Eq. (6.42). 

First-order momentum continuity equation for baryons 

At first order we find 



drj 



+ W1)'+<I>(1)'^ = -T 



4 p^ 



3 pb 



„(!)« _ ^;(1)^ 



(6.43) 
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Second-order momentum continuity equation for baryons 

At second-order there are various simplifications. In particular notice that the term on the right - 
Eq. (6.42) which is proportional to 5b vanishes when matched to expansion of the left-hand side by 
first-order equation (6.45). Thus at the end we find a very simple equation 



hand side of 
virtue of the 



2 \ OT] or] 






3 Ph 



^;(2)» 



,(2)i^ 



r^K 



2 2 

with r' = ^fie'^TO'- 

First-order momentum continuity equation for CDM 

Since CDM particles are collisionless, at first order we find 



(6.44) 






d-q 



Hv^ 



CDM 



(6.45) 



Second-order momentum continuity equation for CDM 



At second-order we find 



1 / dv, 
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dr\ 
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2^ + 2Hc^, + <i>(2)'M- 
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"CDM + "CDM '-'j "CDM + 1^ +* )^ + 



TcDM 
m-CDM 



= 0. 



(6.46) 



VII. ENERGY MOMENTUM TENSORS 

In this section we provide the expressions for the energy-momentum tensors for photons and massive particles in 
terms of their distribution fucntions. 



A. Energy momentum tensor for photons 

The energy momentum tensor for photons is defined as 






cPP Pf'P,, 



■f, 



~^J (27r)3 PO 
where g is the determinant of the metric (2.1) and / is the distribution function. We thus obtain 
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(7.2) 
(7.3) 
(7.4) 
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where p~f is the background energy density of photons and 

are the quadrupolc moments of the photons. The photon velocity is defined in Eqs. (6.31), (6.32) and (6.33). 

B. Energy momentum tensor for massive particles 

The energy momentum tensor for massive particles of mass m, number density n and degrees of freedom g^ is 



-9 
where g is the distribution function. For electrons (or protons) gd = 2 (we are not counting antiparticles) . We obtain 

T^ - -P™ - -Prn [l + <5,« + i(5(2)^ , (7.7) 

T^ = -e^'^'^Prnvl = -e*+*p-„ (^.(,1)^ + i.L^)^ + 5(^).(^)^) , (7.8) 

^m j = Pm i^'j-^ + ^In^^ jj = Prn i^' J -^ + ^'m^^^'m^' j • (7-9) 

where pm is the background energy density of the massive particles and we have included the equilibrium temperature 
Tm- The velocity are defined in terms of the distribution function in Eq. (6.13). 

VIII. SUMMARY 

In this paper we took the first step towards the evaluation of the full radiation transfer at second order in perturbation 
theory by computating the second-order Boltzmann equations describing the evolution of the baryon-photon fluid. 
They allow to follow the time evolution of the CMB anisotropics at second-order and at all scales from the early epoch, 
when the cosmological perturbations were generated, to the present through the recombination era. The dynamics 
at the second order is particularly important when dealing with the issue of NG in the CMB anisotropics: many 
mechanisms for the generation of the primordial inhomogeneities predict a level of NG in the curvature perturbation 
which might be detectable by present and future experiments. Having an accurate theoretical prediction of the CMB 
anisotropy NG in terms of the primordial non-Gaussian seeds is somewhat mandatory. This paper will be followed 
by a companion one where we will present the computation of the three-point correlation function at recombination 
making use of the set of Boltzmann equations derived here. 
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Symbol 



Definition 



Equation 



UJi 

Xij 
V 

f 

9 

f(i) 

Ai) 
Jim 

C(P) 

p 

Ue 
5e{St) 

k 

''CDM 

Sim 
T 



Gravitational potentials in Poisson gauge (2-1) 

2nd-order vector perturbation in Poisson gauge (2-1) 

2nd-order tensor perturbation in Poisson gauge (2-1) 

Conformal time (2-1) 

Photon distribution function (3-8) 

Distribution function for massive particles (4-3) & (6.2) 

i-th order perturbation of the photon distribution function (3.23) 

Moments of the photon distribution function (4.17) 

Collision term (4.2) & (4.6) 

Magnitude of photon momentum (p — pn^) (3.3) 

Propagation direction (3.6) 

First-order fractional energy photon fluctuations (5.5) 

Second-order fractional energy photon fluctuations (5.10) 

Electron number density (6.12) 

Electron (baryon) density perturbation (4.13) 

Wavenumber (5.17) 

Baryon velocity perturbation (5.21) & (5.22) 

Cold dark matter velocity (6.46) 

Second-order photon velocity (6.33) 

Temperature source term (5.18) 

Optical depth (5.15) 

Background photon (baryons) energy density (6.45) 
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